When a positively charged particle passes close to an atom, its path will be deflected (or scatter) by an angle $\theta$ which obeys the relation:
$$\tan (\theta/2) = \frac{Z e^2}{2 \pi \epsilon_0 E b} $$where $Z$ is the atomic number, $e$ is the electric charge, $\epsilon_0$ is the permittivity of free space, $E$ is the kinetic energy of the incident particle, and $b$ is the impact parameter (see diagram).
This process is called "Rutherford Scattering" after Ernest Rutherford, who was among the first physicists to explore this process.
We will model a beam of 1 million $\alpha$ particles with a 2D Gaussian profile incident on a single atom of Gold. We would like to calculate the fraction of these particles which "bounce back" (ie. scatter through angles greater than $\theta = 90^\circ$.
When bounce back happens, therefore $\tan (\theta/2) \gt 1$, so
$$ b \lt \frac{Z e^2}{2 \pi \epsilon_0 E} $$Write a programme which simulates the incident Gaussian beam and calculates the fraction of particles which bounce back.
Please write your own function to calculate Gaussian random numbers using this format:
def gaussian():
YOUR ALGORTHIM HERE
return x,y
Make use of these parameters for the incident beam:
Hint: you can make use of the astropy.constants module to import various constants you need for this problem.
In [19]:
from astropy import constants as const
import numpy as np
import matplotlib.pyplot as plt
#This just needed for the Notebook to show plots inline.
%matplotlib inline
In [2]:
print(const.e.value)
print(const.e)
In [27]:
#Atomic Number of Gold
Z = 72
e = const.e.value
E = 7.7e6*e
eps0 = const.eps0.value
sigma = const.a0.value/100.
#print(Z,e,E,eps0,sigma)
N = 1000000 #Start small, and increase to 1 million when you're sure the code runs correctly.
#Function to generate two sets of random Gaussian numbers.
def gaussian():
r = np.sqrt(-2*sigma*sigma*np.log(1-np.random.random()))
theta=2*np.pi*np.random.random()
x=r*np.cos(theta)
y=r*np.sin(theta)
return x,y
#Main Programme
count = 0 #Initate count of particles bounced back
for i in range(N):
x,y=gaussian()
b=np.sqrt(x*x+y*y)
#If this is true the particle is bounced back
if b<Z*e*e/(2*np.pi*eps0*E):
count +=1
print(count, "particles were reflected out of ", N, "incident")
print("this is a bounce fraction of {0:.5f} +/- {1:.5f}".format(count/N,np.sqrt(count)/N))
Notice something about $b$?
In [28]:
#Atomic Number of Gold
Z = 79
e = const.e.value
E = 7.7e6*e
eps0 = const.eps0.value
sigma = const.a0.value/100.
#print(Z,e,E,eps0,sigma)
N = 1000000 #Start small, and increase to 1 million when you're sure the code runs correctly.
#Main Programme
count = 0 #Initate count of particles bounced back
for i in range(N):
b= np.sqrt(-2*sigma*sigma*np.log(1-np.random.random()))
#If this is true the particle is bounced back
if b<Z*e*e/(2*np.pi*eps0*E):
count +=1
print(count, "particles were reflected out of ", N, "incident")
print("this is a bounce fraction of {0:.5f} +/- {1:.5f}".format(count/N,np.sqrt(count)/N))
In [29]:
?np.random.normal
In [33]:
#Atomic Number of Gold
Z = 79
e = const.e.value
E = 7.7e6*e
eps0 = const.eps0.value
sigma = const.a0.value/100.
print(Z,e,E,eps0,sigma)
N = 1000 #Start small, and increase to 1 million when you're sure the code runs correctly.
#Main Programme
count = 0 #Initate count of particles bounced back
for i in range(N):
x=np.random.normal(0,sigma,1)
y=np.random.normal(0,sigma,1)
b=np.sqrt(x*x+y*y)
#If this is true the particle is bounced back
if b<Z*e*e/(2*np.pi*eps0*E):
count +=1
print(count, "particles were reflected out of ", N, "incident")
print("this is a bounce fraction of {0:.5f} +/- {1:.5f}".format(count/N,np.sqrt(count)/N))
Write a programme following the method of Monte Carlo Integration to calculate
$$ I = \int_0^2 \sin^2 [\frac{1}{x(2-x)}] dx. $$As you will need to calculate $f(x) = \sin^2 [\frac{1}{x(2-x)}]$ many times please write a user defined function for this part of your programme.
In [32]:
#Define the function
def f(x):
fx = (np.sin(1/(x*(2-x))))**2
return fx
#Integrate the function from x=0-2
#Note that you need to know the maximum value of the function
#over this range (which is y=1), and therefore the area of the box
#from which we draw random number is A=2.
N=1000000
k=0
for i in range(N):
x=2*np.random.random()
y=np.random.random()
if y<f(x):
k+=1
A=2.
I=A*k/N
print("The integral is equal to I = ",I)